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Abstract 

Reanalysis of Einstein IPC data and new observations from the GINGA LAC indicate the 
presence of extended X-ray emission (10-50 kpc) around the starburst galaxy M82. Here we 
model this emission by calculating numerical hydrodynamic simulations of the starburst event 
to much later times and larger scales than previously considered. For our models, we adopt a 
supernova rate of 0.1 yr~^, and an extended low-density static halo that is bound to the galaxy. 
There are three stages to the evolution of the wind-blown bubble and the propagation of the 
shock front: the bubble expands in an almost uniform density disk gas, with a deceleration of 
the shock front (t< 3.6 Myr); breakout from the disk and the upward acceleration of the shock 
front (3.6 Myr <t< 18 Myr); propagation into the halo, leading to a more spherical system 
and shock deceleration (18 Myr <t). For a halo density of 10~^cm~^, the outflow reaches a 
distance of 40-50 kpc from the center of the starburst galaxy in 50 Myr. We calculate the time 
evolution of the X-ray luminosity and find that the extended starburst emits 3 x lO^^erg s~^ 
to 10^°erg s-i in the GINGA LAC band and ~ lO^^erg s'^ in the Einstein or ROSAT HRI 
band. The degree of the ionization equilibrium in the outflow and its effect on the iron Ka 
line emission are discussed. 

Key words: Starburst galaxies — Supernova explosions — X-ray emissions — Hot gas 
halos 
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1. Introduction 

Spiral galaxies with enhanced rates of star formation in their central regions have come 
under increasing scrutiny in an effort to understand the causes and evolution of this activity. 
These systems, often referred to as starburst galaxies, possess considerable masses of molecular 
and atomic gas and are among the brightest sources of far infrared emission, presumably 
reradiated ultraviolet light from many early-type stars. Based on the infrared luminosity, the 
CO mass, and the presence of many radio continuum point sources, it is estimated that a 
starburst galaxy such as M82 (NGC 3034) is producing massive stars that become supernovae 
at a rate of ~ 0.1 SN yr~^ (Kronberg and Sramek 1985). With this high rate of supernovae in 
a relatively small volume, the supernova remnants will merge quickly, creating a large bubble 
of hot gas in the center of the galaxy. 

X-ray observations of starburst galaxies such as M82 and NGC253 reveal the presence of 
strong diffuse X-ray emission. The Einstein Observatory HRI observations reveal that the 
emission extends perpendicular to the disk 2-3 kpc for M82 (Watson, Stanger, and Griffiths 
1984; Kronberg, Bierman, and Schwab 1985) and 1-1.5 kpc for NGC 253 (Fabbiano and 
Trinchieri 1984). The extended nature of this emission is indicative that the hot gas has been 
able to flow away from the central region of the galaxy. Hydrodynamic models were developed 
to understand this outflow process in more detail and to reproduce the X-ray observations 
(Chevalier and Clegg 1985; Tomisaka and Ikcuchi 1988). For various adopted initial gas dis- 
tributions, Tomisaka and Ikeuchi (1988) calculated the evolution of the hot gas component, 
where in the early stages, the hot medium is confined to a bubble. Eventually, this bubble 
elongates perpendicular to the disk (along the steepest density gradient) and begins to flow 
outward in a columnated fashion; in these calculations, the hot gas is followed to a distance of 
2-3 kpc from the midplane of the disk. These results appear to be successful in that Einstein 
Observatory X-ray data can be understood within the context of the model. 

Recent observations regarding the extent of the X-ray emission have provided surprising 
results that indicate our understanding of the starburst phenomenon is incomplete. Einstein 
IPC observations (Fabbiano 1988) indicate that fainter X-ray emission extends to about ~ 10 
kpc from the center of M82, but recent observation with the Ginga X-ray satellite revealed 
the unexpected result that weak diffuse X-ray emission extends to distances of several tens of 
kiloparsecs (Tsuru et al. 1990). It is this extremely extended X-ray emission, which was not 
anticipated in the published models, that we investigate here with a series of gas-dynamical 
models that follow the starburst phenomenon to late times and include such effects as an 
extended halo (or intergalactic) medium. Although these calculations are aimed at under- 
standing the recent Ginga observations, they are presented in a more general fashion so as to 
be applicable to a broader range of issues. 

In the next section, the model of the galaxy and the initial density distribution are presented 
as well as the adopted numerical method. In section 3, numerical results are described, while 
section 4 is devoted to the discussion of the expected X-rays emission, including the properties 
of the iron Ka line (/li/ ~6.7 keV). 
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2. The Model and the Numerical Method 



2.1. The Initial Density Distribution 



As in Tomisaka and Ikeuchi (1988; hereafter Paper I), we assume that the gravitational 
field is represented by a spherical stellar component with a modified King density distribution 



[l + (^)2]3/2' 



(2.1) 



where w is the distance from the center of the galaxy. The core radius and the stellar 
density at the center pc are determined to fit the rotation law. The gravitational potential is 
expressed as 



ln^- + [l + (-)T/' 



(2.2) 



We assume that the gas distribution possesses cylindrical symmetry and use a cylindrical 
coordinate system [z, r, (p) to describe the hydrodynamics of the gas flow. For the initial 
conditions, we assume that a part of the radial component of gravity e'^{d(f)^:/dr) (e < 1) is 
balanced by a centrifugal force v^/r, with the remainder of the gravitational force balanced 
by a pressure gradient. If one were to simply assume that the gas rotates rigidly on cylinders, 
then the resulting hydrostatic isothermal density distribution would be 



p = po exp 



(2.3) 



However, a difficulty with this distribution is that at large distances (tu ^ Vc) the gravity is 
weak, so the balance between the centrifugal force and the inward-directed density pressure 
gradient produces a strong funnel around the 2;-axis (see Paper I) which may not be physically 
realistic. Consequently, we relax the constraint leading to this result by introducing the factor 
exp[—{z/Z(i)'^], thus allowing for slower rotation in the halo than the disk, an efi^ect that is 
observed in M82 by Sofue et al. (1992). Then, the resulting initial density distribution used 

here is - 9 9 ^ 

(t)*{z,r) - e^exp [-{z/zdY\ 0*(O,r) [ 



p = po exp 



(2.4) 



The various parameters that define the gravitational potential are chosen so as to reproduce 
known observations. To fit the CO rotation curve of Nakai (1986), we use the parameters 



Tc ^ 350pc, 



(2.5) 
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M.^4w?-10»M,(^)^(3g^)e-. (2.6) 

To include the presence of a halo or intergalactic gas, we assume the presence of a hydrostatic 
hot isothermal component of the form 

/ 0.(^,r) -e^exp [-{z/zaf] 0.(0, r) | 
Ph = Pm exp < 2 r ' (2-7) 

where Cg^ and p^Q are the sound speed and central gas density of this hot component; the 
total density distribution is shown in figure 1. On the z-axis, the cool (disk) component 
dominates below ~ 2 kpc, above which the hot (halo) component is distributed. Here we take 
Cs — 30km s~^, Csh — 300km s~^, po = 20cm~^, — 2x 10~^cm~^, Zd — ^ kpc and e = 0.9. 
With this high value of Csh relative to the gravitational potential, ph is almost constant. Due 
to the high halo temperature ~ 6.5 x lO^K and the low density ~ 2 x 10^^cm~'^, the cooling 
time in the halo - 2 x 10V(^/6.5 x 10^K)(n/2 x lQ-^cm-^)-^{K/2 x lO^^^erg cm^ s~^)~^, is 
longer than the typical time scale of the computation <10^ yr. 



2.2. Basic Equations and Numerical Method 

The two-dimensional, cylindrical symmetric, hydrodynamical equations are the basic equa- 
tions, which are given by 

d pv \ Op 

+ g-^{pv.v,) + -g-^irpv.vr) = + pg:, (2.9) 

d pv d \ d dp 



with 



P 



(2.12) 
7 - 1 

where A represents the radiative cooling coefficient (we use the cooling function derived by 
Raymond, Cox, and Smith 1976), g* is the effective potential that includes the effect of the 
centrifugal force, and the other symbols have their usual meanings. 

As the rotation of the gas was shown to be unimportant (Paper I), we do not solve the 
equation of motion for the rotational component Vy,. The effect of the rotational velocity is 
to reduce the gravity, leading to the effective gravity g* = Vp/p given for the initial gas 
distribution. The r-component of this effective gravity coincides with the ordinary effective 
potential as 

s; = - i . (2.13) 

or r 
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although it should be noted that g* is no longer irrotational. 

The "monotonic scheme" (van Leer 1977; Norman and Winkler 1986) was adopted to solve 
the hydrodynamical equations, i.e., equations (2.8)-(2.12), where the advection and source 
terms are calculated separately. For numerical stability at the shock front, a tensor artificial 
viscosity is included in equations (2.9)-(2.11) (see Tomisaka (1992)). A coohng function is 
fitted to the Raymond-Cox-Smith function using a third-order spline interpolation. To include 
the radiative cooling, the source step of the energy equation (2.11) is solved by two steps, 
where the first step is 



At 

and the second step is 

- e 1 



= —work done by artificial viscosity, (2-14) 



At 2 



(pW + p("+i))divv(") - n(") 2^(T("+i)), (2.15) 



where the superscript refers to the time-step; the cooling function is rewritten as A(T, n) — 
n^L(T). Since the (n + 1) quantities enter the second step in a nonlinear fashion (implicit 
scheme), this equation must be solved iteratively to obtain the new thermal energy e*^'^"'"^^ = 
pin+i) — fikT^'^'^^^ /(7 — 1). Although we employ an implicit scheme to solve the energy 
equation, this scheme is totally vectorized and runs very efficiently on supercomputers. 



3. Results 



3.1. Characteristics of Flow 

The numerical results of outflow arc calculated using the conditions believed to prevail in the 
central region of M82, rsN = O.lyr^^, and uq = 20cm^'^ (other parameters are given in section 
2). The time evolution of the gas density, pressure, and velocity are shown in contour maps in 
figure 2, while the quantitative distribution of these quantities along the 2;-axis are presented 
in figure 3. As the calculation progressed, a larger scale was required, so at early times (figure 
2a) the cell size was Az = Ar = 25 pc and the number of cells was (n^.n,.) = (500,250), at 
intermediate times (figures 2b and 2c) the cell size was Az = Ar = 50 pc and grid size was 
{uz, rir) = (500, 250), while at later times, Az = Ar = 50 pc and the number of cells was 1000 
by 500. 

From the previous calculations on the outfiow from starburst nuclei (Paper I) and super- 
bubbles driven by supernova explosions (for review see Tenorio-Tagle and Bohdenheimer 1988 
and Tomisaka 1991), a hot spherical bubble is initially produced at the center of the starburst 
but soon it begins to elongate in the z-direction when it becomes larger than 2-3 times the 
cold gas scale- height. The effect is clearly seen in figure 2a and in Paper I (figures 4a, and 4b); 
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the cold gas half thickness is related to Zd- During this stage of the evolution, the usual double- 
shock structure is found (figure 3a) where matter ejected from stellar winds and supernovae 
pass through an inward- facing shock (Weaver et al. 1977). The swept-up material pushes into 
the surrounding region, driving an outward-going shock into the cool material. This shock is 
radiative so the hot bubble is surrounded by a cooled shell of material. 

The structure of the hot gas begins to change at an age of about t = 10 Myr, when the 
upper part of the bubble encounters the low-density halo region (n ~ 2 x lO^^cm"^ and 
nT ~ 10"^cm~^ K), where the radiative cooling time is much longer than the expansion time of 
the bubble. Furthermore, the shock begins to accelerate upward when it enters the low density 
material, and this leads to fragmentation of the cool shell, as is seen in figure 2b, where the 
open shell has a funnel shape. Because of the long cooling time, a new cool shell at the top of 
the flow does not form, although there is still a cool shell in the lower part of the bubble that 
is in the disk. At this stage, the upward flow is well-columnated and the gas moves almost 
parallel to the z-axis at supersonic speeds (flgure 2b). The flow is supersonic along the z-axis 
for z>j 1 kpc, it decelerates at the interface between the bubble and the undisturbed medium, 
and then flows outward, which contributes to the tangential growth of the bubble. The other 
contributor to the tangential expansion is the non-negligible azimuthal velocity component of 
the flow from the disk. 

The essential condition that leads to a bipolar outflowing structure is the initial density 
stratification perpendicular to the galactic disk. This was shown in Paper I by comparing the 
model with an initial density distribution similar to that adopted in the present paper to a 
model with an exact plane-stratified density distribution. Although our initial conditions help 
promote a bipolar structure, the bipolar fiow is nearly as pronounced in the calculations with 
a plane-stratified distribution. 

In the next phase of the starburst phenomenon (18 Myr, figures 2c, 3c), the structure and 
dynamics of the hot gas are dominated by the extended low-density medium. A new feature, 
which is common in simulations of radio jets, is the appearance of shocks along the z-axis 
that rcprcssurizc and columnate the fiow. At the base of the flow, from {z,r) ~(Okpc, Okpc) 
to (5kpc, Okpc), the fiow experiences free supersonic expansion, but at a height of 5kpc, the 
pressure of this fiow region falls below that of its surroundings and a shock ensues, reheating 
the gas to >10'^K. Not only does a Mach disk appear, but a shock also occurs along the sides of 
the fiow and may be traced upward from [z, r) ~(Okpc, Ikpc) to the Mach disk near 2; ~ 5 kpc. 
This shock is refiectcd at the symmetry z-axis via Mach refiection and the weak refiected shock 
can be traced further upward. Following this first Mach disk, the fiow expands once again, 
losing pressure and leading to another (weaker) repressurization shock along the z-axis. This 
second shock can be traced from {z,r) ~(10kpc, Okpc) to (5kpc, Skpc), with a sharp change 
of direction near (lOkpc, 2kpc). At heights beyond lOkpc, the jet moves at constant velocity 
and with little expansion in the azimuthal direction. This part of the fiow {z = 10-20kpc) is 
supersonic with respect to the hot halo. As a result, since the gas fiows upward supersonically, 
another outward-facing shock is formed in the ambient intergalactic matter and inside the 
shock front the ambient matter accumulates (see figure 3c; the distribution of mach number 
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for Vz clearly indicates several shocks and it shows that the flow at -2>10 kpc is supersonic). 
There is also a signiflcant azimuthal component to part of the flow, which is caused by the 
flow through the oblique shocks that come off the Mach disks. This component of the flow is 
partly responsible for azimuthal expansion of the hot bubble, making it more nearly spherical. 

In the final stage of development, at 50 Myr (figures 2d and 3d), the disturbance triggered 
by supernova explosions reaches a height of z ~ 40 kpc from the starburst galaxy. At this 
time, the hot gas would entirely envelope the disk of the galaxy. As at the preceding time, 
two repressurizing shocks occur, although further out from the origin, and much of the flow is 
along the z-axis, but less narrowly focused. 

These simulations possess three dynamical phases of the hot bubble that are associated 
with its evolution from the disk into the extended low-density halo. By examining the height 
of the top of the outer-most shock front, Zs, against the time elapsed after the supernova 
burst began (figure 4), we see that the shock front decelerates during the first phase when 
the expansion occurs in the uniform disk material (d'^Zg/dt'^ < for t^3.6Myr). As the shock 
front pushes out of the disk and encounters a steep density gradient and low density material 
perpendicular to the disk (3.6Myr<t<18Myr), it accelerates and it is at this time that the 
initial dense shell fragments. The final deceleration phase occurs because the bubble begins 
to expand in the halo that has a uniform density (t^lSMyr). The final expansion rate of this 
enormous hot bubble is similar to that of a pressure-driven superbubbles and is by given by 
Zs oc t^^^ (Weaver et al. 1976). If the fiow were spherically symmetric, the size would be 

i?~31kpcQxlO-W3) (3xl042erg s'l) (sOM^) ' ^^"^^ 

where L, and Uhaio mean, respectively, the mechanical luminosity of supernova explosions 
E^sn'^sn, and the number density of the intergalactic medium. However, the fiow is actually 
non-spherical and is focused towards the ^-axis. The height of the actual outfiow is larger 
than that expected assuming spherical symmetry {R). Nevertheless, the geometric mean of 
the height of the top {zs) and the half width (rmax) of the bubble {zsTraax)^^"^ agrees well with 
the above estimate [R) and it shows that the expansion is proportional to oc t^/^. Although the 
outfiow shows the bipolarity due to the initial gas distribution, the dynamics of this process 
seems similar to the wind-blown bubble. Previous calculations (e.g., Paper I) have studied the 
bubble development as far as the acceleration phase, but the late stages of development have 
not previously been investigated in the starburst context. 

3.2. The Effect of Halo Pressure and Density 

To understand the effect of the physical state in the halo to the outflow, models with a 

range of different halo parameters were calculated. Model B has a halo with a temperature 
and pressure that is ten times smaller than in model A but with the same density distribution. 
The gravitational fleld is balanced by the pressure gradient. 



7 



Figure 5 shows the structure of this model at t —50 Myr. Compared with figure 2d, we 
see that the overall structure is similar in the two models. The expansion in the ^-direction 
is almost identical in the two models, a result that is probably due to the importance of the 
upward-moving momentum and that the halo density distributions are the same. However, 
the expansion in the r-direction shows differences in that the maximum radial size in model 
B is ~ 0.8 times smaller than in model A. Comparing figure 5b (cross-cut view of figure 5a) 
with figure 2d, the pressure near the expanding front of this model is approximately a factor 
30% smaller than that of the previous model, which is a result of differences in the initial halo 
temperatures in the two models. Thus, due to the decrease in the driving force, the expansion 
in the r-direction is slower than in the previous model. 

To investigate this further, we calculated another model (C), in which the halo has the same 
density distribution of previous models but has a temperature and pressure that is 10 times 
higher than in model A. Once again, the expansion in the z direction is largley unaltered, with 
model C showing a shghtly faster expansion than model A (10%) at t =20My. In contrast, the 
radial expansion is much faster than model A. At t =20Myr, the maximum radial size of the 
outfiow reaches ~20 kpc, while in model A it is only ~12 kpc. 

In model D we have investigated the effect of changing the halo density by calculating a 
model with the same pressure as in model A but with one-tenth the density. 

As expected, the outfiow in a low-density halo expands more quickly than model A once the 
bubble has expanded into the halo (t>5Myr: figTirc 4). While the expansion in the 2;-dircction 
for models A and B are almost identical, the deceleration in the final phase becomes inefficient 
due to the lower density in model D. 

Prom these models A-D, we concluded that the expansion in the 2;-direction is sensitive to 
the halo density but less sensitive to the halo pressure. In contrast, the expansion in the r- 
direction is affected by both the halo pressure and density in the sense that the radial expansion 
becomes faster with higher halo pressure or with lower halo density. 

4. Discussion 

4-1- Expected X-rays — Comparison with Observation 

In this subsection, we calculate the expected X-ray emission by the hot gas outflowing from 
starburst galaxies. We assume collisional equilibrium in the outflow and use the emissivity 
calculated by Raymond et al. (1976). Assuming that the X-ray emission is optically thin, the 
total luminosities of various X-ray bands are calculated by 



where Aeand is a band X-ray emissivity coefficient while ^out and Tout are the sizes of the 
numerical outer boundaries. In figure 6, the expected X-ray luminosity for model B is shown 




(4.1) 
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as a function of time for three energy band: 0.155-24.8 keV (solid line), 0.284-1.56 keV 
(dotted line), and 1.56-8.265 keV (dashed line). It is assumed that the 0.284-1.56 keV band 
approximates the Einstein Observatory HRI band and that of 1.56-8.265 keV approximates 
the Ginga LAC band. Fore- and background emission were subtracted from the value derived 
by equation(4.1), assuming that the expected X-ray luminosities at the initial state should 
vanish. (Before the supernova explosions begin, the expected luminosities in the HRI and 
LAC bands from the galaxy and the extended halo are estimated to be 1.7 x 10^^ erg s~^ and 
6.8 X lO^^erg s~^, respectively.) 

Before 5 Myr, the X-ray luminosities decrease, an effect that is due to the progressively 
decreasing density as the bubble expands and the shock front moves upward (this phase was 
studied in detail in Paper I). In the subsequent ~5 Myr, the X-ray luminosities become a 
factor of 5 - 10 brighter than the minimum value at t ~ 5 Myr. This brightening occurs when 
the shock begins to accelerate {t ~ 3.6Myr) which leads to a higher post-shock temperature 
and a higher luminosity. After 10 Myr, the LAC band X-ray luminosity stays nearly constant, 
while in the lower energy bands, the luminosity increases slowly. The fact that, at late times, 
the X-ray luminosities do not decrease monotonically but remain fairly constant is important 
observationally in that one expects to see this phenomenon over a long period of time, although 
at progressively lower surface brightnesses. 

To more deeply understand the evolution of the X-ray luminosity, we first note that the 
X-ray emission is a combination of the shock-heated matter originally ejected from the parent 
galaxy and the halo material compressed by the outward-facing shock front. If we adopt a 
simple spherically symmetric approximation for the evolution of the material, this hot gas 
drives a shock into the halo gas. The time averaged temperature of the shock-heated material 
ejected from supernovae is Tj„ = {l/5)V^{fimp/k), where denotes the wind velocity and it 
is equal to = {2rsNEsN / MsnY^"^ (Weaver et al. 1977). This temperature, which remains 
constant during the course of the starburst event, is Tin — lO^K for parameters of our models 
A and B. (This temperature seems to depend upon the geometry: in the actual starburst, a 
one-dimensional outflow like a funnel is formed and the temperature is below lO^K, while in 
the spherical symmetric configuration the temperature must be ^ lO'^K. This is because in an 
oblique shock the fraction of kinetic energy in the post-shock region is larger than that of the 
perpendicular shock.) The shock speed decreases with time as it progresses further into the halo 
gas, so the post-shock temperature of the compressed halo matter decreases monotonically. In 
addition to the evolution of the hot starburst event, the observed X-ray emissivity also depends 
on the temperature behavior of the bandpass cooling function (Raymond et al. 1976): Ahri 
is a decreasing function of T for T>10^K and remains constant for 2 x 10®K<T<10^K, while 
Alac is an increasing function of T for T<10^K. Most of the emission of the HRI band seems 
come from the compressed halo matter. Since the shock is decerelated and the post-shock 
temperature is decreasing from ~ lO^K (after 10 Myr), the emission in the HRI band increases 
due to the increase of Ahri. On the other hand, since the emissivity decreases after ~ 10 Myr, 
the luminosity of the LAC band does not increase but remains constant, even if the amount 
of the accumulated mass increases. 
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Tsuru et al. (1990) reported that the luminosity emitted from the extended region sur- 
rounding M82 in the 2-10 keV band is approximately ~ 4 x 10^^ erg s~^, which corresponds to 
~ 10 % of the total luminosity from M82. In comparison, the X-ray luminosity calculated here, 
3 X 10^^ — 10"^° erg s~^, is consistent with the observed value. Furthermore, the size of the X-ray 
emission region reaches ~ 50 kpc in half width and the size is expected to continue to grow 
without decreasing X-ray luminosity, provided that there is a halo or intergalactic medium 
with the same density at these larger sizes. This extended X-ray emission is consistent with 
the observational width as seen by the LAC, which is wider than the point-spread function of 
LAC (FWHM of LAC is 1.1 degree and it corresponds to the hnear size of 57(i:)/3.25Mpc) 
kpc, where D represents the distance of the object). Using our calculations, the mass of the 
hot gas is estimated to be Mhot ~ 1.5 x lO'^M0(r/5OMyr). 

The mass loss rate of model B corresponds to the assumption that the amount of ejected 
mass per supernova explosion AM is equal to 30 Mq. This value of the mass includes not only 
the ejecta of supernovae, but other sources, such as mass loss from stars that are not supernova 
progenitors or evaporated mass from clouds engulfed in the hot outflow. We adopted this mass 
loss rate in models A-D. However, we could have taken the limiting case that the only mass 
entering the system was the supernova ejecta, in which case the value of AM would be ~ 5M0. 
We have calculated a model (E) using this lower mass loss rate, in order to understand how the 
X-ray luminosity would be affected; the remaining parameters are identical to those in model 
B. Since decreasing the mass loss rate leads to a higher input gas temperature, the energy flux 
into the system remains constant. If a significant component of the X-ray emission were to 
come from this ejected matter, we might expect the X-ray spectrum to be harder. 

The time variations of X-ray luminosities are shown in figure 7. Llac decreases from 
3 X 10^^ erg s~^ (model B) to 10'^^ erg s~^ (model E), while Lhri remains ~ 5 x 10^*^ erg s~^, 
which seems to contradict the above expectation. The X-ray luminosity in the HRl band 
must come largely from the accumulated halo/intergalactic matter, which depends on the 
total energy input to the system but is insensitive to the mass and temperature of the ejected 
matter. If the X-ray emission in the LAC band is also from the accumulated intergalactic 
matter, the luminosity will be expected to be the same as that of model B. Thus, a part of the 
luminosity in the LAC band is emitted from the ejected matter. Since the ejected mass-loss 
rate is taken to be 1/6 that of model B, the emission measure of the ejected matter seems to 
become smaller than that of model B. Even if we take into account that the emissivity Alac 
becomes higher than that of model B, the total luminosity in the LAC band becomes smaller 
than from model B. 

4-2. Iron Line Emission 



4.2. a. Iron Line Emission from the Starburst Galaxy 
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Unfortunately, due to the low intensity of the extended emission, the spectrum of the 
extended emission around M82 was not detemined by the LAC on Ginga (Tsuru et al. 1990) 
and must be left to the next generation of X-ray telescopes, ASTRO-D, XMM or AXAF. 
However, the spectra of the X-ray emission from the central part of nearby starburst galaxies, 
such as M82(Makishima and Ohashi 1989; Ohashi et al. 1990a) and NGC 253(Ohashi et 
al. 1990b), were obtained by Ginga. These investigators fit optically thin thermal emission 
models to M82 and NGC 253 and obtained temperatures Tx — 5.7 ± 0.3 keV (M82) and 
6.0 ± 0.8 keV (NGC253). They report equivalent widths for the iron emission line at 6.7 keV 
of 170±60 eV(M82) and 180±?^g cV (NGC253). It came as a surprise that there is no definitive 
evidence of iron line emission in either system, because if the X-ray emission is optically thin, 
thermal, and comes from the hot gas created by supernova explosions in starburst galaxies, the 
iron Ka emission line hu ~ 6.7 keV should have been detected. The expected equivalent width 
is approximately ~ 900 eV, if we assume that the coUisional ionization balance is achieved and 
that the gas has a temperature determined by fitting the observed continuum spectra (Ohashi 
et al. 1990a). The failure to detect the iron emission may indicate that the iron abundance 
is less than 1/5 - 1/4 of solar abundance or that some non-equilibrium ionization process has 
invalidated some of the underlying assumptions. 

To estimate the expected Fe abundance, we note that a type II supernova progenitor with 
an intial mass of 20 Mq ejects ~ O.IM© (Thielemann, Hashimoto, and Nomoto 1990) of iron. 
As adopted in our models, this amount of iron becomes mixed with 30 Mq of material, leading 
to a metalicity that is several times the solar value. The amount of material that becomes 
mixed with the ejecta (30 Mq) cannot be significantly larger (300-500 Mq would be needed to 
properly reduce the observed Fe abundance) or it would create an X-ray continuum luminosity 
in excess of that observed. 

We also consider whether including non-equilibrium ionization effects can affect the iron 
abundance inferred from the observations. Although non-equilibrium ionization models have 
not been developed for starburst galaxies, such models have been developed for supernova 
remnants, for which the physics is quite similar (Itoh 1977, 1978, and 1979; ShuU 1982; Hamil- 
ton, Sarazin, and Chevalier 1983). From hot plasma, the X-ray emission at energies hv< 
IkeV comes from the line emission of He- like low-Z ions such as C, N, and O. Higher-energy 
X-ray emission is dominated by the continuum emission from Hydrogen and Helium (free-free 
emission). When ionization equilibrium is not achieved (i.e., the atoms are in an ionization 
state lower than expected under equilibrium conditions), low-Z ions emit above their equi- 
librium value. In contrast, the emissivities of high-Z ions, such as Fe and Ni, become lower 
than their equilibrium values. This difference occurs because, for temperatures near ~ 10^ K, 
the equilibrium ionization state of these low-Z species is higher than for the He-like species. 
The numerical result of Hamilton et al. (1983) show that the luminosity of Fe (and Ni) Ka 
line can be less than 0.1 times that expected from a plasma in ionization equilibrium. The 
departure from the equilibrium values depends upon the ionization time scale nor, where uq 
and T represent the ambient density around the supernova remnant and time elapsed since the 
explosion. For constant values of UoT (see their equation 4b and figure 6), we see that when 
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TT-oT = lO^yr cm~^, the ratio of emission in the non-equihbrium stage to that in the colhsional 
equihbrium stage, F, is approximately equal to 0.3 for a shock temperature of Ts ~ lO^K. In 
an earlier phase, where uqt = 500yr cm~^, F has a value of ^0.1. Thus, it appears that if the 
gas has not had sufficient time to achieve ionization equilibrium, the 6.7 keV iron line emission 
may be much lower than value expected from ionization equilibrium, which might explain the 
small equivalent widths of iron Ka emission lines in starburst galaxies. 

If the average gas density in a starburst nucleus is ~ 10~^cm~^, the ejected matter travels 
a distance ~ Ikpc before achieving ionization equilibrium (t ~ 1 Myr). This situation is 
much different from the interstellar space in nomal galaxies, where all the supernova remnants 
except for extraordinary young remnants with an age of < 1000 yr achieve the ionization 
equilibrium. The ionization equilibrium may have been achieved in a lower halo or a disk 
and the outflowing matter may be in the equilibrium. However, the X-ray emission from the 
starburst galaxy seems to be different from that expected from the emission of supernova 
remnants in our galaxy. 

4.2.b. Iron Line Emission from the Starburst Outflow 

Although the spectrum of the X-ray emission coming from the starburst outflows has not 
been observed, here we estimate the strength of the iron hue of the large-scale outflow. We 
estimate the effective "nor" for the starburst outflow by integrating the quantity 



along the z-axis. The factor oc means that the value should be volume averaged. For the 
Sedov solution, this average value equals 4.17noT. This quantity gives us a mean measure of 
how close the plasma is to a state of equilibrium ionization. Here, since in starburst galaxies 
the ejected gas seems to escape more freely than ordinary disk galaxies, we assume that the 
ejected gas has a lower ionization state than expected from its temperature using colhsional 
equilibrium when it leaves the galaxy. 

We plot equation (4.2) as a function of time for model E in flgure 8. Since the undisturbed 
medium is static, the value ndz/v^ goes to infinity in the pre-shock medium. The gradually 
increasing part corresponds to the outflow material. From this flgure, we see that the integrated 
value is large ~ 1.2 x 10^ yr cm~^ during the early phase {t —5 Myr). But after the outflow 
breaks out of the gas disk completely (>8 Myr), the swept up material is of much lower 
density, so the integrated value decreases and is not more than ^ 6000 yr cm~^. Since the 
ejected matter after the break-out can flow with slight interaction of ambient material, the 
integrated value decreases ~ 6000 yr cm~^. This low value seems to indicate a low degree of 
equilibriation. In a late phase, the value becomes large again due to its large size and long 
age. 

If we use the conversion factor from the value of equation (4.2) to rioT for the Sedov solution, 
the value ~ 6000yr cm~^ corresponds to ~ 1500 yr cm~^ for uoT. This implies that the material 




(4.2) 
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outflowing from a starburst galaxy shows modest non-equilibrium effects that could affect the 
interpretation of the iron abundance based upon the Fe Ka line. If we estimate the ratio of 
the emission in the non-equilibrium to the equilibrium state (F) based upon this value of nor 
as applied to the supernova remnant calculations (T ~ lO^K) of Hamilton et al. (1983), we 
find that F = 0.1 — 0.3. More exact simulations, i.e., 2-D hydrodynamical simulation with the 
non-equilibrium ionization, are necessary to more accurately answer the problem, a calculation 
we hope to undertake in the near future. 
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ment, University of Michigan, Ann Arbor, for their hospitality and kindness, where he began 
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National Astronomical Observatory, and Facom VP200 at Institute of Space and Aeronautics 
Sciences. This work was supported in part by the Grant-in-Aid (02854016) from the Ministry 
of Education, Science and Culture and by grant NAGW-2135 from National Aeronautics and 
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Tables 



TABLE 1: Model Parameters 



Model 




MsN 


rico 


riho 






Csh 


spacing 


number of cells 






(lo^^g s-^-: 


(cm-'') 




) 


(km s-i) 


(km s-"^) 


(pO 




A ... 


0.1 


1.9 


20 


2 X 10" 


-3 


30 


300 


50x50 


500x250 
1000x500 


Al .. 


0.1 


1.9 


20 


2 X 10- 


-3 


30 


300 


25x25 


500x250 


B ... 


0.1 


1.9 


20 


2 X 10" 


-3 


30 


95 


50x50 


1000x500 


C ... 


0.1 


1.9 


20 


2 X 10" 


-3 


30 


950 


50x50 


500x250 


D ... 


0.1 


1.9 


20 


2 X 10" 


-4 


30 


950 


50x50 


500x250 


E ... 


0.1 


0.3 


20 


2 X 10" 


-3 


30 


95 


50x50 


500x250 
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Figure Captions 

Fig.l: (a) The initial state of the density (left panel) and pressure (right panel) distributions. 
The disk component occupies the region below za, while an extensive halo of nearly contant 
density lies above it. (b) The initial distribution of pressure (dash-dotted line) and density 
(dashed line) are plotted along the 2;-axis (left panel) and the r-axis (right panel), where 
the standard pressure, p*, is taken 10^ K cm~^. 

Fig.2: The structure of the starburst outflow for model A at four times: t — 5Myr (a), 
t — lOMyr (h),t — 20Myr (c), and t — 50Myr (d). The right panel shows pressure contours 
and velocity vectors while the left panel shows density contours. Before t =5Myr, the bubble 
is already focused by the effect of the initial density gradient, and after t =8Myr, the bubble 
expands in a uniform- density halo, producing shocks in the process. 

Fig.3: Cross-sectional views of the starburst outflow along the z-axis (left panel) and the 
r-axis (right panel), showing velocity log|f^/f*| (solid line), temperature logT(K) (dash- 
dotted line), density logn(cm^'^) (dashed line), and Mach number = Vz/{'yp/p). Here, 
= 11.65km s~^, is the isothermal sound speed for lO^K gas. Figures (a) - (d) correspond 
to (a) - (d) of figure 2, respectively. The heights of the cross-section views (right panel) are 
a,t z — Okpc (a), 5kpc (b, c), and lOkpc(d). 

Fig. 4: The expansion law of the outflow, showing the time variation of the outermost shock 
front traveling on the 2;-axis. This clearly shows that the expansion is divided into three 
phases. First, a hot bubble is created in the high density disk, and the shock decelerates. 
Second, when the size of the bubble surpasses the size of the disk, it elongates in the z- 
direction and the shock acclerates into the stratifled density distribution. Finally, the shock 
begins to decelerate when it encounters the uniform-density halo. 

Fig. 5: (a) The structure of the starburst outflow for model B at t = 50Myr, where the 
right panel shows pressure contours and velocity vectors while the left panel shows density 
contours. This model has the same density distribution as model A but has a pressure 
and a temperature 10 times smaller. Comparing this with flgure 2d, the expansion in the 
2;-direction is almost identical, but the expansion in the r-dircction differs by ~ 40%. (b) 
Cross-sectional views of the starburst outflow along the z-axis (left panel) and the line 
z — lOkpc (right panel), showing velocity log\vz/v:^\ (sohd line), temperature logT(K) 
(dash-dotted line), density logn(cm~^) (dashed line), and Mach number = Vz/i'jp/ p)- 

Fig. 6: The time variations of X-ray luminosity from the outflowing hot gas (model B) are 
given in the energy bands 0.155-24.8 keV (solid line), 0.284-1.56 keV (dotted line), and 
1.56-8.265 keV (dashed line). The second and third band correspond, respectively, to those 
of the Einstein Observatory HRI and the Ginga LAC. 

Fig. 7: The same as flgure 6 but for model E. Here, the mass loss rate from the supernovae 
plus evolving stars is 1/6 that of model B; all other parameters are the same as in model B. 

Fig.8: The distribution of Jq ndt along the 2;-axis for model E. Since the initial halo material 
is static, the value of ndt is artificially large at early times. However, after 8 Myr, the 
value of ndt is never more than ~ 6000 yr cm~^, so that the emitting gas is in a modestly 
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non-equilibrium state. This may help to explain the low observed value for the equivalent 
widths of the iron Ka emission lines. 
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